Understanding critical collapse of a scalar field 



Carsten Gundlach 

LAEFF-INTA (Laboratorio de Astrofisica Espacial y Fisica Fundamental - Instituto Nacional de Tecnologia Aerospacial) , 

PO Box 50727, 28080 Madrid, Spam 
(8 April, 1996) 



I construct a spherically symmetric solution for a massless real scalar field minimally coupled to 
general relativity which is discretely self-similar (DSS) and regular. This solution coincides with 
the intermediate attractor found by Choptuik in critical gravitational collapse. The echoing period 
is A = 3.4453 ± 0.0005. The solution is continued to the future self-similarity horizon, which is 
also the future light cone of a naked singularity. The scalar field and metric are C 1 but not C 2 
at this Cauchy horizon. The curvature is finite nevertheless, and the horizon carries regular null 
data. These are very nearly flat. The solution has exactly one growing perturbation mode, thus 
\Q ' confirming the standard explanation for universality. The growth of this mode corresponds to a 

Q\ ' critical exponent of 7 = 0.374 ±0.001, in agreement with the best experimental value. I predict that 

in critical collapse dominated by a DSS critical solution, the scaling of the black hole mass shows 
a periodic wiggle, which like 7 is universal. My results carry over to the free complex scalar field. 
$H Connections with previous investigations of self-similar scalar field solutions are discussed, as well 

as an interpretation of A and 7 as anomalous dimensions. 
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I. INTRODUCTION 



■ A. Critical phenomena in gravitational collapse 

In an astrophysical context, gravitational collapse normally starts from a star. This means that the initial data are 
almost stationary, and that they have a characteristic scale which is provided by the matter. Therefore astrophysical 
black holes have a minimum mass, namely the Chandrasekhar mass. Abandoning the restriction to almost stationary 
initial data, or alternatively to realistic matter, one should be able to make arbitrarily small black holes. One may 
then ask what happens if one tries to make a black hole of infinitesimal mass by fine-tuning the initial data. 

The investigation is simplified by choosing a matter model that does not admit stable stationary solutions. Then, 
for any initial data, there are only two possible outcomes: formation of a black hole, or dispersion leaving behind 
flat spacetime. The first systematic numerical examination of the limit between the two (the "critical surface" in 
phase space) was carried out by Choptuik [Q for a massless minimally coupled real scalar field in spherical symmetry. 
He evolved members of various one-parameter families of initial data each of which comprised both collapsing and 
dispersing data, and searched for the critical parameter value by bisection. For all families he investigated he was 
able to make arbitrarily small black holes by tuning the parameter p of the data: there was no evidence for a "mass 
gap" . Instead he found two unexpected new phenomena. 

For marginal data, both supercritical and subcritical, the time evolution approaches a certain universal solution 
which is the same for all the one-parameter families of data. This solution is an "intermediate attractor" in the sense 
that the time evolution first converges onto it, but then diverges from it eventually, to either form a black hole or 
to disperse. This universal solution (also called the "critical solution") has a curious symmetry: it is periodic in the 
logarithm of spacetime scale, with a period of A ~ 3.44. This is also referred to as "echoing" , or discrete self-similarity 
(DSS). 

Moreover, for marginally supercritical data, the final black hole mass scales as M ~ (p — p*) 7 , where p is the 
parameter of the family of initial data, and its critical value. The "critical exponent" 7 has the value ~ 0.37 for the 
scalar field, and like the critical solution it is universal in the sense that it is the same for all one-parameter families 
of data. 

Both phenomena were then also found in the axisymmetric collapse of pure gravity Q , indicating that they are an 
artifact of neither the choice of matter nor of spherical symmetry. There, A was found to be ~ 0.6, and the critical 
exponent 7 to be ~ 0.36. For a perfect fluid with equation of state p = p/3 in spherical symmetry ||, the universal 
attractor has a different symmetry: it is not discretely, but continuously self-similar (CSS). 7 is found to be ~ 0.36 
once more. 

Choptuik's results have been duplicated, to varying precision, in ^-|). 
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Subsequently, the matter models were generalized. For a fluid with p — kp in spherical symmetry, now with 
arbitrary constant k, 7 was found to be strongly fc-dependent []. The real scalar field model was generalized 
to a one-parameter family of two-component non-linear sigma models l|. This family includes the cases of a free 
complex scalar field ||[l0|], a real scalar field coupled to Brans-Dicke gravity JTTJ and, as a special case of the latter, a 
string-inspired axion-dilaton model fl2|| . For the axion-dilaton model, 7 ~ 0.264 is found in collapse simulations Jl3| . 
From these new examples it is clear that 7 is not universal with respect to different kinds of matter, but only with 
respect to the initial data for any one matter model. 

B. The emerging picture 

Let us now examine some general features of the "critical solutions" which appear as intermediate attractors in 
collapse, and which seem to describe the limiting case of the formation of a zero mass black hole. First of all, they 
must be scale-invariant in some way, and in fact will show homotheticity, or "self-similarity of the first kind" Jl4) . 
Because of self-similarity, they must have a curvature singularity, but they should not have an event horizon. The 
absence of a horizon means that the collapsing matter is visible in the solution. 

The unique endpoint of gravitational collapse is given by the Kerr-Newman family of solutions because, roughly 
speaking, they only admit damped perturbations: they are attractors in phase space. If the critical collapse solution 
were also attractors, we would see many naked singularities in nature, and this is clearly not the case. In fact they 
are attractors of co-dimension one, and we shall see that this gives rise to an analogue of the "no hair" theorem: 
universality. The Kerr-Newman and critical collapse solutions are briefly contrasted in Table |[ 

In a schematic picture of phase space Jl5[ , the critical solution remains within the critical surface. The observed 
universality suggests that it is in fact an attractor within the critical surface, and in consequence an attractor of 
co-dimension one in phase space. This attractor could either be a fixed point, corresponding to CSS, or a limit cycle, 
corresponding to DSS. Fig. [l] illustrates this. 

Nearly critical Cauchy data are situated close to the critical surface, but may be far from the critical point. Under 
time evolution they are attracted towards the critical point. While they approach it, their "distance" from the critical 
surface increases exponentially but remains small until they are close to the critical point because, by the assumption 
of near-criticality, it is initially very small. Near the critical point, the exponential increase takes over, and the phase 
space trajectories are repelled from the critical surface, all into the same direction. This constitutes the mechanism 
of universality with respect to initial data. The formula for the black hole mass follows essentially by dimensional 
analysis, and the critical exponent 7 can be read off from the linear perturbations of the critical solution |15|. As we 



shall see in section [VB ( the periodicity of the critical solution in the DSS case gives rise to a periodic wiggle in the 
scaling law. 

What confuses the naive phase space picture is gauge-invariance: the same spacetime corresponds to very different 
trajectories in superspace, depending on how it is sliced. The naive picture must therefore be used with care pending 
its formulation in geometric terms. 

CSS solutions and their linear perturbations have already been calculated for various matter models. For the 
p = p/3 model ||,[l5| and for the axion-dilaton model Jl2|,[l3| they agree with the experimentally determined critical 
solution and give the correct critical exponent. For the complex scalar field a CSS solution has been found ||, but it 
is only an attractor of co-dimension three For the p = kp model with k ^ 1/3, 7 has been calculated formally 
on on the basis of the perturbations of a CSS solution but it is not known for what range of k this CSS solution 
is really the critical solution. (By continuity, it must be for some neighborhood of k = 1/3.) Furthermore, no CSS 
solution exists for k > 0.888 M, and one would expect that the critical solution becomes DSS at either that, or a 
smaller, value of k. The CSS solution and its perturbations have been calculated also for the family of nonlinear 
sigma models M. The number of unstable modes of the CSS solution changes from one to three at some value of 
the parameter. This probably indicates the changeover from CSS to DSS in the critical solution. Why some critical 
solutions are CSS and others DSS is not yet understood, however. 

The present paper gives the first calculation of a DSS critical solution, together with its maximal extension and 
its linear perturbations. This is technically much more difficult than CSS, but DSS is the most generic case of self- 
similarity, and the mathematical and numerical methods developed here will be useful in other applications. The 
critical mass scaling generalizes to contain a universal periodic wiggle. 



1 The 7 for k 7^ 1/3 have been calculated only in linear perturbation analysis, but I strongly expect them to be confirmed by 
collapse calculations. 
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The plan of the paper is as follows. In section O, I define DSS and construct the DSS solution of the real scalar 
field model in the past light cone of the singularity as a nonlinear eigenvalue problem. It agrees with the critical 
solution found by Choptuik B| . This section, together with appendices [a|, |c[ ^ and [f| is an expansion of the Letter 



L6[. Sections Jni] and |V| both build on section ||, but are independent of each other. 
In section 



[II, I extend the critical solution up to the future light cone of the naked singularity. I find that this 



Cauchy horizon is regular, and that it is in fact nearly flat. In section [V, I calculate the linear perturbations, and 
generalize the calculation of the critical exponent to the DSS case. My value for 7 agrees with the experimental one, 
but I also predict the existence of a (small) universal wiggle overlaid on the power-law scaling, which has not been 
observed so far. I show that the Choptuik solution is the critical solution not only for the real, but also the free 
complex scalar field. In section |v|, I put my results into the context of results for other collapsing systems on one 
hand, and of the study of (continuously) self-similar scalar field models on the other. I then discuss the next steps 
to be taken, and the possible connection with critical phenomena in statistical mechanics and quantum field theory. 
Various details are given in the appendices. 

II. THE CRITICAL SOLUTION 
A. Field equations 

In this section I construct an isolated solution of general relativity minimally coupled to a massless real scalar field 
with the following properties: 1) spherical symmetry, 2) discrete self-similarity (to be defined below), 3) analytic at 
the center of spherical symmetry, 4) analytic at the past self-similarity horizon, 5) the scalar field is bounded. (It is 
possible that there is no solution obeying 1) to 4) that does not also obey 5), but I have not shown this.) 

The Einstein equations we consider here are 

Gab = &kG (^<j),a(j),b - -^gab4>,c(t>' C ^j , (1) 

in spherical symmetry. The matter equation <f> c ' c — follows from the Einstein equations as the contracted Bianchi 
identity. The spacetime metric is 

ds 2 ee -a(r, t) 2 dt 2 + a{r, i) 2 dr 2 + r 2 (d9 2 + sin 2 9 dip 2 ). (2) 

This form of the metric is invariant under transformations t — > t(t), a — > a, such that adt = adt. In order to write 
the matter equations in first order form, we introduce the auxiliary matter fields 



X(r, t) ee sph^G -<j> r , Y(r, t) = V^rG -<f> t . (3) 
a ' a ' 

The combinations X± — X ± Y of these fields propagate along characteristics. The radial null geodesies, which are 
also the matter characteristics, are characterized by the quantity j = o/a |. The scalar wave equation in X± is then 

r (X ± . r t gX±. t ) = [±r^ - r^] X ± -X T . (4) 
la a J 

In the following we use X + , X_ , g and a as our basic variables. A complete set of Einstein equations in these variables 



is 



2 

and the matter equations become 



r 9.r= (1 - a 2 )g, (5) 
ra r = -a [(1 - a 2 ) + a 2 {X\ + X 2 )] ee C u (6) 

gra^]-a z (X 2 + ~ X 2 _) ee C 2 , (7) 



2 I have made one change of notation. In jlfj I denned g = exp(^o)a/a, while g = a/a here. This is for greater convenience in 
the remainder of the paper. 
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r (X± lT T gX± t , 



-(l-a 2 )-a 2 X 



X± - X T = c± 



(8) 



when we eliminate the metric derivatives with the help of the Einstein equations. The five first order equations (|5||8|) 
are our field equations. I have defined the expressions C±, C\ and C2 for later convenience. The absence of g it in the 
equations reflects the fact that a, and hence g, contains a gauge degree of freedom not determined by the Cauchy 
data. 

The two scalars made from the Ricci curvature, using the Einstein equations, are 



R = 4r- 2 X+X. 



r> nao 
Kabtt 



R 2 . 



(9) 



The Riemann tensor will be considered in appendix 



B. Discrete self-similarity 



The concept of (continuous) self-similarity (CSS) (or homotheticity) has been defined in a relativistic context JL7 14 
as the presence of a vector field x such that 

£ x 9ab = 2g a b, (10) 

where C x denotes the Lie derivative. I now introduce the concept of discrete self-similarity (DSS). In this symmetry 
there exist a diffeomorphism c/> and a real constant A such that, for any integer n, 

(</>*)" 9ab = e 2nA g ab , (11) 

where c6* is the pull-back of <f>. 

To see what DSS looks like in coordinate terms, we introduce coordinates (a, x a ), such that if a point p has 
coordinates (cr, x a ), its image (f){p) has coordinates (cr 4- A, a;"). One can verify that DSS in these coordinates is 
equivalent to 

g^(a,x a ) = e 2CT .g MI/ (cr,x Q ), where g^{a, x a ) = g^(a + A, x a ) (12) 

In other words, the DSS acts as a discrete isomorphism on the rescaled metric g^v a is intuitively speaking the 
logarithm of spacetime scale. 

One can formally construct such a coordinate system in the following way: Fix a hypersurface £ such that its image 
£' under <p does not intersect E. Introduce coordinates x a in S, and copy them to £' with <j). Introduce coordinates 
(cr, x a ) in the region between £ and £' such that a has range [0, A], their restriction to £ is (0, x a ) and their restriction 
to £' is (A, x a ). Finally, copy these coordinates to the entire spacetime, such that if p has coordinates (cr, x a ), its n-th 
image c6"(p) is assigned coordinates (cr + nA,x a ). Clearly there is enormous freedom in defining such coordinates. 

In order to clarify the connection between CSS and DSS, one may define a vector field x = d/da, although there is 
no unique x associated with a given <j). The discrete diffeomorphism 4> is then realized as the Lie dragging along x by 
a distance A. Clearly, CSS corresponds to DSS for infinitesimally small A, and hence for all A, and is in this sense a 
degenerate case of DSS. In this limit, x becomes unique. In the following I speak of DSS only in the absence of CSS. 

In order to see what form DSS takes in spherical symmetry in the particular coordinates defined in (||) , we make a 
coordinate transformation t = e a T(a 1 z) and r = e c ' R(u, z), where T and R are periodic in cr withperiod A. Here cr 
is the same as in the general construction, and x a = (z, 6, ip). In the new coordinates the metric (0) takes the form 

ds 2 = e 2a j-a 2 [(T + T a ) da + T z dz] 2 + a 2 [(R + R <a ) da + R. z dz} 2 + R 2 (d9 2 + sin 2 9 dtp 2 ) j (13) 

This is of the form ( p"2| ) if and only if a(z, a) and a(z, a) are also periodic in a with period A. In terms of t and r 
this periodicity corresponds to 

a(r, t) = a{e nA r, e nA t), a(r, t) = a(e nA r, e nA t). (14) 

This is not yet the most general way to impose DSS in (|^). We obtain that by also admitting the transformations 
t — > i(t), with a — * a — dt/di a. a need no longer be periodic, but it must be related in this way to some a that is. 
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C. Formulation as an eigenvalue problem 



We now introduce specific coordinates of the kind just discussed. The following choice will be sufficiently general 
for our purpose: 



In" ' 



C = ln(I)-£ (r). (15) 



Here r is an arbitrary fixed scale, and £o( T ) is a periodic function with period A^|. Both are to be determined later. 
In the new coordinates, the matter and Einstein equations are 

X±X= l±(l + Qe<^ 9 (16) 
g x = (1 - a 2 )g (17) 

o,C= °i ( 18 ) 
o,r= e-K+f^ff" 1 ^ + (1 + Qd (19) 

Here £ = d£,o{r)/dT, and in this paper a prime always denotes the derivative of a function of one variable with 
respect to its formal argument. The equations are invariant under a translation in r, corresponding to a change in 
the arbitrary scale r . 



In order to impose discrete homotheticity, we demand the periodicity (14). We also impose the regularity condition 
a = 1 at r = and the gauge condition a = 1 at r = 0. (Both are compatible with the periodicity.) In our choice of 
dependent and independent variables we therefore impose boundary conditions 

a(C,T + nA) =a(C,r), g((, r + nA) = a(C, r), (20) 

and 

o(C = -oo,r) = l, <?(C = -oo,r) = 1. (21) 

Note that we continue to describe the metric with the variables a and a, but that they are not the coefficients of the 
metric associated with coordinates (r, Q . 

From the Einstein equations it follows that the periodicity condition must hold also for X + and X_. From the 
equations defining X and Y, we obtain 

lT = (2irG)~ 1/2 a \(1 + QX + g^e-^+^y] , jC = (27rG)~ 1/2 aX (22) 

Because the right-hand sides of both equations are periodic in r, the scalar field itself is of the form 

4>((, t + A) = (periodic in r) + kt, (23) 



where k is a constant, k is not an independent parameter, but is determined through the first of equations (|22j). <p is 
bounded and periodic if and only if k = 0. 

As we have a 1 + 1 dimensional hyperbolic problem, we can interchange space and time. In this view, near r = 0, 
equations ([l6]), ( p"7| ) and (|lj) form a first-order system of evolution equations for a, g, X + and X_ with "time" 
coordinate C and (periodic) "space" coordinate r. The data are subject to one constraint (|19|), which is propagated 
by the evolution equations. As £ — > — oo, we impose boundary conditions corresponding to a regular center r = 0. 

The equations become singular when the denominator of (Hq) vanishes. The treatment of this singular surface, both 
analytical and numerical, is simplified if we use the coordinate freedom incorporated in the choice of £o( T ) to make 
this happen "for all r at once", namely on the line C — 0- We therefore impose the coordinate condition 

[l-(l+£ )eS] f=0 = 0. (24) 



3 As in pq |, I assume t > 0. In a collapse context, where the spacetime region we are about to calculate is to the past of the 
singularity, t then decreases to the future, but this is purely a matter of convention. The convention can be reversed in the 
results simply by changing the sign of Y. 
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A necessary condition for regularity is now that the denominator also vanishes at £ = 0, or 



[C_ - e> 9 X_, T ] = 0. 



(25) 



We may look at this from a different angle: The coordinate condition means in fact that ( = is null, and hence a 
characteristic of the equations. Such a characteristic which is mapped onto itself by the self-similarity map <f> is called 
a "self-similarity horizon" (SSH), or a "sonic line" |l8). For our equations, with periodic boundary conditions in the 
"space" coordinate r, it constitutes a Cauchy horizon. On every other £ = const, surface, X + and X_ would be free 
Cauchy data, but on £ = they are overdetermined, as one requires only characteristic data. This gives again rise to 
the condition (p5|). 

We still have to see what kind of regularity this condition enforces, and if it is sufficient as well as necessary. To 
investigate the possible behavior of the solution at the past SSH £ = 0, we assume for the moment that a, g and X + 
are at least C° there, and expand the equation for X_ to leading order in The resulting approximate equation, 
using the coordinate condition (E4T) and equation (p7|), is 



A(t)X_ - X_ tT + C(t) 
0B(t) 



where the coefficients 



(1 - a 2 Q ) - a 2 X_ 



D 



are evaluated at C = 0. This approximate equation admits an exact general solution, namely 

x = x» lhom ( T ) + X^ om (C, r). 
The particular inhomogeneous solution X 1 ^ " 1 is defined as the unique solution of 

Ax inhom _ ^inhom + C = Q 



(26) 



(27) 



(28) 



(29) 



with periodic boundary conditions. This solution exists and is unique, unless the average value of A vanishes. The 
general homogeneous solution X^° m is of the form 



X 



horn 



C 



4a i a-£p-i b 



J o e 



F 



IqB - In |C| 



(30) 



The notation here is that of appendix Aq is the constant part, or average value, of the periodic function A. IqA 
denotes the principal function of A — Aq, that is of A minus its average, with the integration constant defined so 
that IqA itself has zero average. F is a periodic function of one variable with period A. It is the free parameter of 
the solution. We do not need to determine it here. What is important is that the solution is analytic if and only if 
F vanishes identically. If F does not vanish, there are two possibilities. If Aq/Bq < 0, the solution will blow up at 
C = 0. If Aq/Bq > 0, the solution will be Cq but not C\ there. It is easy to see that both Aq and Bq are negative in 
our case, so that the solution either blows up or is analytic. So with the one condition ( p5| ) we automatically enforce 
analyticity. 

We shall impose one extra symmetry on our ansatz. The results of numerical collapse simulations jlj] indicate that 
<f> itself is periodic. Moreover, even if one adds a potential to the scalar field action, the attractor is found to be the 
same as for the massless field JO^. This requires that 4> remains bounded in the critical solution, because only then 
can a polynomial of 4> be neglected with respect to the space and time derivatives of 0, which are unbounded because 
of the echoing on an exponentially decreasing spacetime scale. For this reason we look for a solution with k = 0. 
Because all frequencies are coupled through the nonlincarity of the field equations, K vanishes if and only if the even 
frequencies of X± and the odd frequencies of a and g vanish for all £. Therefore we now impose in our ansatz that 
a and <?, as well as £o> are composed only of even frequency terms in r, and X± only of odd frequency terms. This 
reflection- type symmetry is compatible with the field equations. 

We have now completed the formulation of a hyperbolic boundary value problem. Its Cauchy data are values for 
the four fields a, <?, X + and up to one constraint, and up to a translation invariance in t, plus the unknown value 
of A and the unknown function £o( T )- The count is therefore 4oo — oo — 1 + 1 + oo = 4oo. (Here oo stands for the 
countably infinite number of degrees of freedom of a periodic function of one variable.) These free data are balanced 



G 



by two boundary conditions at £ = — oo and two at £ = 0, or 4oo degrees of freedom again. One would therefore 
expect this boundary value problem to have at most a discrete set of solutions. 

Numerically I have constructed one such solution. Locally it is unique. To solve the boundary value problem 
numerically, I have expanded all periodic fields in Fourier components of r, truncating the expansion at some relatively 
small number of components. This takes advantage of the fact that the solution is smooth, so that the Fourier series 
converges rapidly. The field a is not evolved in £, but reconstructed at each step from the constraint Details 
are given in the appendices. I find that A = 3.4453 ± 0.0005. In a previous publication fTql , I had given a value of 
A = 3.4439±0.0004. The difference, corresponding to 3.5 standard errors, is due to a change in numerical details of the 
algorithm. These changes correct what I now regard as an inconsistency in my original Pseudo-Fourier method, and 
are therefore "systematic error" . The quoted error is in both cases estimated discretisation error, which is discussed 
in appendix [F]. The present Fourier methods are outlined in appendix 

I have shown in that the DSS solution I have constructed agrees with the intermediate attractor observed by 
Choptuik jlj] to the numerical precision of the latter. The error in the DSS solution I have corrected here is small 
enough not to affect this agreement. 



III. MAXIMAL EXTENSION OF THE CRITICAL SPACETIME 



A. From the past self-similarity horizon to t = 

The coordinates ((,t) become singular at t = (£ = oo). Clearly it will be necessary to replace r with p as the 
periodic coordinate a to regularize the equations there. Before we do this, it is useful to examine the asymptotic 
behavior of the solution in the old variables as £ — * oo, or t — > 0. Neglecting terms of order exp(— £), the field 
equations in this limit reduce to 

(l+QX±, c =X ±tT , (31) 
(1+Qa c =a r , (32) 
a c = Ci, (33) 

g,c= (i - « 2 ).9, (34) 



with the general solution 



X±= X ±0 (p), (35) 
a= a (p), (36) 
g=g (r)e^ + ^, (37) 



where 



p = C + T + £o(T)=ln(^j , (38) 

and the periodic functions X±q and ga are free parameters of the general solution. The periodic function ao is uniquely 
determined by X± through the ODE 

ao = ~ao[(l-a§) + ag(X£o + X* )], (39) 
and the the periodic function a and constant v are determined from ao by integration, as 

i/ = (l^§) , <T = J„(l-a§), (40) 

where the notation is that of appendix 

From this asymptotic expansion we see that a and X± are regular at t = 0, while g is not and will have to be 
replaced by another dependent variable. 

We replace £ and r by p as above and w given by 

w = e ( 1 + n )(i—p)+f(T)+ h (p)^ (41) 
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where n is a constant and h(p) and /(r) are periodic functions of one variable, all yet to be determined. Among the 
dependent variables, we replace g by 

g = [1 + n + /'(r)] e «(^-p)+/M+Mp) ff . (42) 

We now have DSS if the dependent variables a, g, .X+ and X_ are periodic in p at constant w, with the same period 
A we have determined previously. 

The field equations in the new dependent and independent variables are 

X± , w = C± ~ X± f , (43) 



g. p + (a 2 - 2 + T)g 



(44) 



1 U> 

a,™= ff" 1 ^, (45) 
ap^Ci+r^- 1 ^, (46) 

where 

T(p) = l + n-h'(p). (47) 



We are dealing again with four evolution equations (now in w instead of C) ( |43H45| ) and one constraint (|46|). Singular 
points arise where the denominator of ([l3|) or of (Q) vanishes. 

We first examine the singular point w = 0, corresponding to t = 0. We assume that a is at least C° there, as 
suggested by the asymptotics, and takes the value ao(p). If we set a = ao(/?) m equation (|44|), this approximate 
equation admits an exact general solution, which we can write as 



j(w, P ) = exp <^ - ——h(p)-I Q (a* -2 + j)(p) + F 

1 1 + n 



In \w\ — h(p) 

1+71 



(48) 



Here v is defined by equation (^), and i 7, is a periodic function of one variable with period A, which serves as the free 
parameter of the solution. We see that g will either blow up or vanish at w = 0, unless we impose n — v. Furthermore 
g will not be C° unless F = 0. But imposing these two conditions is equivalent to the one condition 

[g. p + (a 2 -2 + T)g] w ^ = 0. (49) 

(This is so because the expression a 2 — 2 + T can only be the derivative of the periodic function \ng if its average 
value vanishes. This corresponds precisely to n = v, although v does of course not appear explicitly in the equations.) 
This is just the regularity condition one would have expected from inspecting the equations, but now we have shown 
that imposing it actually corresponds to imposing analyticity. 

Before we discuss the other singular points, we use h(p) to identify £ = £q with w = 1, thus simplifying the matching 
between the two coordinate systems. (We then need numerical interpolation only in one, not two dimensions). We 
define the periodic function £o implicitly in terms of £o and Co by the equation 

|o[t + Co+£o(t)]=£o(t). (50) 

(Numerically this is solved by iteration.) Next we define 

/(p) = /[p-Co-|o(p)]- (51) 
This definition implies that / (p) and / (t) coincide when restricted to the line C = Co- Now we fix h as 

h(p) = (1 + n) [Co + lo(p)] - f(p), (52) 

and it can be verified that with this definition w — 1 for C = Co- Now let ao(r) = a(^o,r), and define fio from clq in 
the same way as / from /. Proceed similarly for <?, X + and X_. clq and X±o now constitute initial data for a and 
X± on w = 1 . The initial data for g on w = 1 can be expressed in terms of go as 

S( W = 1, p) = — ^- e Co+ « o(rt ffo(p)- (53) 
1 Co 
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Now we come back to the other singular points of the equations, namely where — Tw _Tg = 0. Straightforward 
calculation shows that this happens where dw = is null, just as the denominator of ( |l6[) vanishes where dC = 
is null. Either occurrence indicates a SSH. If we choose Co large enough, we do not encounter the past SSH. Why 
don't we simply choose Co = 0, however? As one expects, the denominator of the matter equation then vanishes on 
the entire line w — 1, but it vanishes also at two more lines crossing w = 1. This is confirmed by the fact that the 
partial derivative of the denominator with respect to w vanishes at two points (per period) on w = 1. This behavior 
would give us the numerical problems we avoided in the old coordinate system by making the SSH a line of constant 
£. (In contrast, the denominator of the matter equation in the (C, r) coordinate system is everywhere increasing with 
C, which means that it vanishes only at C = and nowhere else.) 

We have now formulated the analytic continuation at t = as a boundary value problem mathematically quite 
similar to the one we have solved above. We consider a as determined by g and X± via the constraint (ff6|). The 
unknowns are the fields g and X± between w = 1 and w = 0, and the function T(p). The three boundary conditions 
on the left are the matching of g and X± to the data at C = Co, and the one boundary condition on the right is (|4S|). 
Note that the unknown functions and constant appear in the field equations only in the one combination T, and that 
g(w = 1) is determined from T, Co and g a . This breaks up our numerical procedure naturally into two steps. In a first 
step, the functions Co, oo, <7o and X±q are once and for all determined from the data at C_= Co- Then we vary T in a 
relaxation algorithm, until we have solved the system (f43|-[46|) with boundary condition (fl9"|). 

We have seen that in the boundary value problem we are only dealing with one free function T(p) to be determined, 
while /, h and n do not appear explicitly. I have introduced them because they are useful in deriving an initial guess 
for r for use in the numerical algorithm. Substituting the asymptotic form (|31]) of g into the definition of g we obtain 

g ~ 5o (t) [1 + n+ /'(t)] e/M^M+^-^^-pH^ri+^P) (54) 

(Note that, exceptionally, I have mixed coordinate systems, in using p and r as independent variables.) This is regular 
at t = oo if n — and if 

5o (r) [1 + n + /'(r)] e /H-««°H = 1. (55) 

This latter condition can be considered as a linear ODE for exp/(r), given go and n — v. By evolving in coordinates 
(C, t) to large C, and comparing with the asymptotic form ( p7|) we can estimate v and go- From there we can calculate 
an estimate for / via (|5^) and then for h and hence T via ( |50| , pl| , |52| ). 

Now that the problem has been solved in the new variables, we can calculate a = a/g from (|f2|). It contains a 
singular i-dependent factor that we can absorb into the (singular) redefinition t — ► t = t 1+n e^ lnt '^ > . The resulting 
regular a is 

g 

n and h are given in terms of the now determined function T by n = T — 1 and h = —IqT. Note that a is no longer 
periodic in p, although the spacetime is DSS. We have had to use the most general form of the metric compatible 
with (|^), abandoning the gauge condition a = 1 at r = 0. In fact a is singular at r — 0, but then we only use it for a 
patch around spacelike infinity. 

We obtain a clearer picture of the behavior at spacelike infinity by ignoring the "wiggles" in a and a. Then they 
simplify to the expressions 

a ~ vl — n, a ~ r~", (57) 

which are valid for r» \t\. The Hawking mass is proportional to the radius, and the geometry at spacelike infinity is 
conical. Our constant n coincides with "1/n" in the notation of Q. Numerically I find n ~ —0.16. 

B. From t = to the future self-similarity horizon 

It remains to extend the spacetime all the way to the future SSH, the future light cone of the singularity at 
(t = 0, r — 0), beyond which the continuation is no longer unique. A priori, the future SSH might itself be singular, 
because we have no free parameters left to make it regular. 

At t = 0, the periodic data X±(p) are fixed, and they determine a(p) through the constraint ([flf). g on the other 
hand is pure gauge at t = 0, depending on our choice of T through equation (f49|). We now continue to evolve in w with 
the same equations as before, but with a new choice of T. This means that we introduce a third coordinate system, 
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although one of the same class ( |3q , pT| ) as before. This time we choose T so as to make the line w = — 1 the future 
SSH, which is also the future light cone of the singularity. Our reason for this is twofold: on one hand we are better 
placed to control the vanishing of the denominator in the matter equation (now for X + instead of A_ ) when it takes 
places on a coordinate line, on the other hand we want the edge of the domain of dependence of our data to coincide 
with a coordinate line, so that we can evolve right up to it. Our boundary conditions are now the data for X± and 
a at w = 0, the constraint (^) which determines g from T at w = (up to a constant factor), and the coordinate 
condition g — T (vanishing of the denominator of A+ jU) ) at w — 1. A priori this new boundary value problem is not 
well-posed. We have no freedom left to impose the vanishing of the numerator of X+ w at w = 1 as an additional 
boundary condition, so that the solution should be genuinely singular at this point. 

To investigate what happens we once more make an analytic approximation. Let us assume that X_ is at least 
C° at w — — 1, and takes value X^. (p). By definition, as our new coordinate condition, g takes the value T(p). In 
consequence X + drops out of the constraint ([l5"|), and independently of the value of X+, a takes the value ao(p), 
which is the solution of the ODE 

ao = 5<«o(l -a§)+<*o*-o (58) 

with periodic boundary conditions. The solution exists and is unique. We can now calculate g >w at w — — 1, and 
hence the linear approximation to the denominator —Tw — g near w = — 1. With these expressions in hand we write 
out the equation for X + , in the leading terms in both the numerator and the denominator. We obtain an approximate 
equation of the same form as (p6|), namely 

_ A{p)X + -X + , p + C{ P ) 

A +,«> — / . 7T7T7 \ 1 l 0y J 

(w + l)B(p) 

where the coefficients are 

A(p) = i(l - 4) - a 2 Xi , B(p) = a 2 -2 + (lnr)', C(p) = -X- . (60) 

The exact general solution to the approximate equation is, once more, 

X + =X i f™\p)+Xl° a \w,p), (61) 

where 

4lf om - X™ h p om + C = (62) 



and 



= + 1)3? e IoA ~t hB F 



I Q B - In \w 
p+- 



B 



(63) 



As we have no freedom to impose any boundary conditions, F does not vanish and the solution is not analytic at 
w = —1. But we see that Bq < 0, while Aq < unless (ao — 1) and X-q vanish identically. The infinitely oscillating 
term therefore vanishes at w = — 1 as (w + 1) £ , where e is positive and small. X + is therefore C°, although not C 1 . 
This is a remarkable result: The presence of even a small amount of radiation crossing the future SSH (the component 
X_) regularizes the radiation running along the horizon (the component X + ), by damping its oscillations. A similar 
result was found by Home in the axion-dilaton [ fl2"| and free complex scalar H (see note added in print) CSS solutions, 
using a similar analytic approximation. 

From the field equations it follows that a, g and X- are C 1 but not C 2 (with respect to w, differentiation with 
respect to p is not a problem). As not all second derivatives of the metric exist, one must ask if the spacetime 
curvature is finite. In appendix ^ I show that all components of the Riemann tensor are in fact C° . 

Although the numerical problem is ill-defined because of the presence of an infinite number of oscillations in X + , 
I have run a naive relaxation algorithm on it. The algorithm does in fact converge, and I even see convergence in 
the values of T, ao and X_ with decreasing step size. I find r ~ 1 + v, and do ~ 1 and A_ — 0. This means that 
spacetime is approximately flat on the future SSH, and very little scalar field radiation is crossing it. X + however, 
oscillates more and more rapidly. These oscillations appear to be at constant amplitude, and I do not see their 
eventual decay numerically. This is consistent with the fact that e -C 1, so that the decay is very slow. 

In order to make the problem numerically well-defined, in spite of the solution being singular, one would have to 
subtract the singular part, solving for F in another boundary value problem. (This caveat may apply also to Home's 
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numerical results I have not attempted this nontrivial step, as I do not see an immediate need for quantitative 

results. Although the problem is not numerically well-defined as it stands, I am confident in the qualitative result 
that the null data at the future SSH are regular, and nearly, but not quite flat. 

Why all this explains in hindsight why the future SSH is, in some sense, regular, it does not explain why it is nearly 
flat, that is, why X_q is so small. I do not see a mechanism in the field equations that would drive arbitrary data 
at w — to very nearly flat space values. And in fact the SSH tends to be less flat when I slightly perturb X + or 
X- at w — (and adjust a accordingly). This indicates that the near- flatness at the future SSH is a property of the 
special DSS solution which is regular at the past SSH. (Incidentally, it is also an argument against near-flatness being 
a numerical artifact.) 



IV. LINEARIZED PERTURBATIONS AND CRITICAL EXPONENT 



A. The linear eigenvalue problem 

Now we turn to the study of linearized perturbations of the critical solution, which leave the perturbed solution 
regular at r — and C — 0. If such a regular perturbation existed while being periodic in r, the critical solution 
would not be isolated. The perturbations must therefore break the discrete homotheticity, that is, the periodicity in 
t. The coefficients of the linearized equations, however, are periodic, and therefore the general linear perturbation is 
a superposition of terms of the form 



8Z((,t) =^C ie A ^Z(C,T), 



(64) 



where each SiZ is periodic in r with period A (although 8Z is not!). Z is a shorthand for (a, g, X + , X-). The 
exponents Xi and hence the 5iZ must be allowed to be complex even for real SZ. As the equations are real, they 
form complex conjugate pairs, corresponding to sine and cosine oscillations in r with frequency QA^. Because the 
ansatz already contains all frequencies that are integer multiples of 2tt/A, and because values of Xi come in complex 
conjugate pairs, we need consider only < 3Ai < 7r/A. 
If the equation for 5Z is of the form 



8Z 



C 



ASZ + BSZ T , 



(65) 



the equation for SiZ is of the form 



5 i Z x = (A + X i B)S i Z + B5 i Z tT . 



(66) 



This indicates how we obtain the equations for SiZ from the linearized equations for SZ. In the following I denote 
the components of SiZ by ^a, Sig and 5iX±. The equations for the periodic quantities SiZ are then 



2 v 2 

a 



S t X ± - (2a 2 A + A_ + 1) S. t X T - aX± (l + 2X\ ) Sitt 



±e C+ «° [X^ T 5 l9 + g (SiX^r + X i S l X±) - (1 + QX ±x S l9 ] \[l±(l + Qe i+ ^g] ' 



fii9X = (1 _ fl2 )^i.9 - 2agSi 



' + la 2 (Xl 
2 2 V + 



Xi - 1 
3 



Sia + a 3 (X + S t X + + X_SiX_) 



S ia , r = -XiSid + e-^+^g- 1 



+ (l+Co){ 



\^{X\ - X*)6ia + a 3 (X + S t X + - X^SiXJ] 



1 3 

2 ' 2 



+ ^a 2 {X\+X 2 _ - 1) 



5,a + a 3 (X+S l X+ + XSiXJ, 



X 2 _)5 t g 



(67) 

(68) 
(69) 

(70) 



We now consider the boundary conditions for these equations. At ( — > — oo we have one free function 8iY\ (r) 
(compare appendices [A| and ^). At the boundary £ = 0, the denominator of the 5iX_ equation vanishes, because it 
is the same as that of the X_ equation. We therefore have to impose the vanishing of the numerator. Imposing the 
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linear constraint as w ell, we can freely specify Sig and 8iX + at £ = 0, and calculate from them Sia and 8iX_. Details 
are given in appendix |B| 

We now have three free functions <J,-Yi, Sigo and SiX-o at the boundaries. At some intermediate point the three 
functions Sjg, 8jX + and SiX- will have to match. The fourth, Sia, will match automatically by virtue of the linearized 
constraint (J70|). We also have the free constant Ai to solve for. Its presence is balanced by the fact that because of 
linearity an overall factor in the perturbations is arbitrary and has to be fixed in some way. The Xi are the eigenvalues 
in a new hyperbolic boundary value problem, this time a linear one. One would therefore expect them to form a 
discrete set. 

Details of the numerical method are given in appendix ^, and an error estimate in appendix [f]. Because of the 
even-odd symmetry of the background, perturbations which have the same symmetry decouple from perturbations 
with the opposite symmetry, and we can consider them separately 

In the left half plane of A^, corresponding to modes that grow towards the singularity r = — oo, I have found one 
perturbation of the first type. In the following I refer to this mode by i = 1. It is real, with Ai = — 2.674 ± 0.009, 
corresponding to a critical exponent 7 = 0.374 ± 0.001. I find no perturbations of the second type in the left half 
plane. There is one perturbation, of the first type, at Ai = 0. It is the gauge mode SZ oc Z*. T , corresponding to a 
translation of the background in r. On the positive real axis I find perturbations of both types, and presumably there 
are more elsewhere in the right half plane. 



B. Critical scaling of the black hole mass 

The derivation of the scaling of the black hole mass for near-critical initial data which I present here was suggested 
in j|, and first made explicit in jl5| for the CSS case. The DSS case requires a subtle generalization, and I find a 
new phenomenon: a "wiggle" in the mass scaling law. 

In preparation, let us consider a family of Cauchy data at constant t, namely 

Z T (r) = Z* ( ln ^' T ) +eSlZ ( ln ^ T ) ' ( 71 ) 

where r is the parameter of the family. Here Z*(£,r) is the critical solution that we have just constructed, and 
e AlT (5iZ(C, t) is the one linear perturbation mode that is growing on small scales, that is, which has negative A. tq is 
the (arbitrary) fixed length scale introduced in (|l5|) and e is a fixed small constant, small enough so that the linear 
approximation is good initially, lnr/ro is the value the formal argument £ takes. Clearly this family is periodic in t 
with period A. 

What happens at late times when these data are evolved in tl For one sign of e, say e < 0, we must find dispersion, 
for the other a black hole. The key observation |l(J is that the data depend on r only through the dimensionless 
combination r/fo, while the evolution equations themselves have no scale. The entire solution scales as 

Z T (r,t)=f(- ) -,r). (72) 

Therefore we know, without having to rely on the linear approximation, that the black hole mass, which has dimension 
length, must be proportional to r$. More precisely, we have 

M = r e" {T \ (73) 

where fJ.(r) is an unknown periodic function of period A. It can be calculated numerically by evolving members of 
the family Z T which span one period in r. 

Now we consider a generic solution. If the initial data are sufficiently close to criticality, there is a spacetime region 
in their evolution where the solution "echos" , that is, where it is close to the critical solution. There the solution is 
well approximated by the critical solution plus linear perturbations, as 

00 

Z(C, t) ~ Z*(C, r) + ]T Ci(p) e^StZiC, r). (74) 

i=l 

The amplitudes Cj of the perturbations depend in a complicated way on the initial data in general and hence on the 
parameter p of a given one-parameter family of initial data. As t — > but while the perturbations are still small, we 
can neglect all perturbations but the growing one, and we have 
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Z(C,r) * Z.(C,t) + d(p) e^SxZiCr). (75) 

By definition we obtain the precisely critical solution for p — p*, and so we must have Ci(p*) = 0. Linearizing C\(p) 
around , we obtain 

Z(C,t)*Z.(C,t) + (p-P*) e^^tCr). (76) 



We define 



1 = ~ (77) 



(as the notation suggests, this will turn out to be the critical exponent), and 

(78) 



t {>>) = - hi ( P - P * ) T] : ; lis 



1 dd . 
(P* 



e dlnp 



Note that T*(p) is only a way of rewriting the "reduced parameter" (p — p*)/p*, while Ti depends on the family of 

initial data through the unknown function C\(p). e is the same as in definition @. If we now fix r = r* (p) + ri in 
the approximate solution (|76|), we obtain a p-dependent family of Cauchy data, namely 

Z P {r) = Z* (in r*(p) + n) + e <SiZ (in r»(p) + , (79) 

where 

r(p) = r e^W+n+foInW+n] ( 80 ) 
But this is of the form for which we have just calculated the black hole mass. Therefore we have 



M(p) = r {p) e ^ ( - p)+Tl] = r (^y^j 



i 

„ti + (m+€o)[t. (p)+Tl] 



(81) 



Let us first consider the CSS case. Then fi and £o degenerate into constants. We recover the well-known exact scaling 
of the black hole mass. The unknown, family-dependent, constant t\ corresponds to an unknown overall factor in the 
mass. 

In the DSS case we find a "wiggle" overlaid on the scaling law, unless the function (fj, + £o) vanishes identically. 
The period of both fj, and £o is nominally A, but £o has only even frequencies, and so does /i, being based only on the 
metric coefficients a and g. The real period in r of (jj, + £o) is therefore A/2, and hence the period of the wiggle in 
the directly measured parameter ln(p — p*) is A/(2j) ~ 4.61. 

The offset of the wiggle is the same constant T\ that already appears in the overall factor. Given the function /x(r), 
the black hole mass is therefore as much determined in the DSS case, namely up to one family-dependent constant, 
as in the CSS case, and \i(t) has the same universal significance as 7. It would therefore be interesting to determine 
/i(r) from evolving ( fn| ) and to test the expression (|8l]) against collapse simulations. 



C. The free complex scalar field model 

As mentioned already, there is a regular spherically symmetric CSS solution for the free complex scalar field, that is 
a complex scalar field with neither a mass term nor a coupling to an electromagnetic field [p| . Later it was discovered 
that this solution has not one but three unstable modes |l(J , and that the critical solution for the free complex scalar 
field is in fact the real DSS solution we have discussed here (up to a global complex phase) p3] . 

The action is, in loose notation, R+ |<9$| 2 . Writing the complex scalar $ as ifi + iip, the action becomes R+ (d(f>) 2 + 
(dip) 2 , or that of two real scalar fields. Any solution of the real scalar field model is therefore also a solution of the free 
complex scalar field model. Furthermore, purely imaginary scalar field perturbations Sip of any purely real solution <j) 
decouple from the purely real perturbations 5<p. This holds because the stress tensor is the some of a term quadratic 
in (j> and one quadratic in ip. The first-order perturbation of this second term vanishes if the background value of ip 
is zero. 

In order to obtain all perturbations of the real scalar DSS solution within the free complex scalar model, I only need 
to calculate its purely imaginary perturbations Sip (in addition to the real perturbations Sep already known). These 
obey the wave equation on the fixed background metric with source , while the corresponding metric perturbations 
vanish to linear order. I have checked that all these modes are damped. With a little extra work I have thus confirmed 
perturbatively that the real DSS solution is an attractor of co-dimension one even in the free complex model. Details 
will be given elsewhere. 
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V. CONCLUSIONS 



A. Critical solutions and matter models 



As I have argued in the introduction, critical solutions of the kind discussed here play a role in critical collapse similar 
to the role the Kerr-Newman solution plays in generic, non-critical, collapse. The crucial difference is the presence of 
matter: A critical solution does not have a horizon, and therefore depends explicitly on the choice of matter. This 
reduces the importance of any one such solution. Two other differences are mainly of technical importance. While the 
Kerr-Newman solutions are known in closed form, all critical solutions discovered so far are only known numerically. 
Furthermore, CSS solutions, like stationary ones, have effectively one dimension fewer, but this simplification does 
not hold for the DSS ansatz: it only makes all fields periodic in one coordinate. The labour involved in constructing 
the real scalar field critical solution, which is DSS, is therefore much greater than for the other self-similar solutions 
found so far, all of which are CSS. There are two reasons for investing it. On one hand the real scalar field model is 
the first in which critical phenomena were observed, and its "echos" , and the "wiggle" in the mass scaling, are new 
features that differ from the well-known CSS case. On the other hand, the real scalar field is a testbed of methods for 
dealing with DSS which can now be applied to other DSS critical solutions. The most interesting of these is probably 
that of pure gravity in axisymmetry |2|] . 

This raises the fundamental problem of angular momentum and electric charge in the initial data for gravitational 
collapse. On one hand, the resulting black hole can have angular momentum and charge. On the other, both must be 
smaller than the mass. So what happens to the black-hole charge and angular momentum if one tries to fine-tune the 
black-hole mass to zero? Clearly the case of rotation is the more interesting one, but it cannot be treated in spherical 
symmetry. Only in one case so far has critical collapse be considered in axisymmetry but unfortunately without 
angular momentum. 

Another restriction of the matter models which have so far been considered in the study of critical phenomena is 
that most of them, in marked difference from any realistic macroscopic matter, do not introduce a preferred length 
scale. (In units where G = c = 1, this is equivalent to the absence of dimensionful parameters in the action.) This 
guarantees the existence of an exactly self-similar solution. Even in the presence of a preferred scale, however, a 
self-similar solution may still be a good approximation towards the singularity, as s 2 = r 2 + t 2 — > 0. This is the case 
for one known exception, namely when a mass-term m 2 <p 2 is added to the scalar field action. Then the DSS solution 
we have constructed here is a good approximation for s <C m , simply because <fi is bounded while <9</> is not. In 
collapse simulations, it is found that the Choptuik solution is in fact the attractor for the massive scalar field |L9| |. 
It remains to be investigated how other matter models introducing a scale react to the attempt to make small black 
holes. A very recent example is the Einstein- Yang-Mills system, which has one intrinsic scale, and which shows both 
a mass gap, and critical mass scaling, for different classes of initial data pi] , A rather different, and interesting, one 
is the attempt to consider semi-classical Einstein equations, with the new scale being the Planck length 



B. Other self-similar scalar field solutions 



In this paper, I have considered a scalar field <p which is (a) bounded on the entire spacetime (k = 0), (b) real, and 
(c) discretely self-similar (DSS). These conditions are suggested by the universal attractor that Choptuik [|J found 
in critical collapse of a real scalar field, and they are justified by the fact that the unique solution of the resulting 
eigenvalue problem coincides with Choptuik's attractor to numerical precision I now review some related work 
on self-similar scalar field solutions departing from one or the other of these assumptions. 

CSS real scalar field solutions have been studied in ^3|. The solution for k = can be derived in closed form and 
was apparently first published in (^4|, then rediscovered in the context of critical phenomena in |2^|2(|. CSS arises 
as a degenerate case of our formalism when a, g and X± are not only pe riod ic in r, but do not depend on r at all. 
In our notation, continuous self-similarity implies that Yi(t) in equation (A2) is a constant, namely Y\ — n. (This is 
incompatible with the even-odd symmetry I have assumed.) 

Brady |^3| has shown that for all k in the CSS case, boundedness at the past CSS is automatic, in marked contrast to 
the DSS case. As in the DSS case, continuation past the SSH is not unique. Brady has considered the one-parameter 
family of possible continuations (for each value of k) . Unfortunately, it is not clear which of his continuations is the 
analytic one. In the DSS case, I have investigated only the analytic continuation at the past SSH, because the past 
SSH is in the Cauchy development of regular data when we think of it as arising as an attractor in collapse simulations, 
and therefore I believe it should be analytic. Of Brady's results I review here only the special case k = 0, because 
it appears to be most closely related to our case. It is qualitatively different from k ^ 0, in that no CSS solution 
with a regular center exists, except for flat empty space. There is a one-parameter family of solutions with a singular 
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center r — 0. One branch of r — is always timelike and has negative mass, the other is either timelike with negative 
mass (called subcritical solutions in J25|), or spacelike with positive mass, and in the latter case its is preceded by 
an apparent horizon (supercritical solutions). In either case both the past and the future SSH carry flat-space null 
data, and one can therefore replace the negative mass part of the solution with flat space in both the past and future 
light cone of the singularity. Subcritical solutions thus pieced together are qualitatively like the maximally extended 
Choptuik solution: r = is a regular center, except at the single point (r = 0,i = 0), which is a naked singularity. 
Both the past and the future SSH are regular and flat. (In the Choptuik solution the past SSH is not flat, and the 
future SSH is only nearly flat.) 

What would happen if we allowed n in DSS, while still imposing regularity at the center and on the past SSH? 
Because of mode coupling, this would mean giving up the even-odd symmetry, and doubling the (infinite) number of 
degrees of freedom in the boundary value problem. Therefore, one might obtain an infinity of new solutions, or no 
new solutions at all. If any new solutions exist, one would wonder next why it is only the one with k = that serves 
as a universal attractor. If there is a family of such solutions continuous with the Choptuik solution, that question 
would be even more acute. I leave these questions to future work. 

A CSS solution (with k = assumed implicitly) has also been constructed for the free complex scalar field [||. In 
this solution, only the metric and the scalar field modulus are CSS, but the complex phase of the scalar field is ilut 
for some constant ui, so that the complex scalar field might be considered DSS with (in our notation) A = 2-k/lo. Not 
surprisingly, our solution shares more features with this one than with the CSS real scalar field models: There is a 
unique solution with a regular center and regular past SSH, and the null data on the future SSH are for nearly flat 
space. The same qualitative picture was found for the CSS solution with axion-dilaton matter |Q. It is remarkable 
that now three different self-similar solutions are known which are nearly, but not quite, flat at the future SSH. The 
Roberts solution is clearly the limiting case, with a flat SSH (but consequently a scalar field that is not C !), and 
may be of help in finding an explanation. 



C. Perturbations, universality, and "renormalisation" 



I have found exactly one unstable mode of the critical solution. The picture of an attractor of co-dimension one is 
thus confirmed perturbatively. The Lyapunov exponent of ~ — 2.674 ± 0.009 gives rise to the value 7 ~ 0.374 ± 0.001 
of the critical exponent. This value agrees with the most precise value obtained from collapse calculations J^], which 
is given as 0.374. Allowing for complex scalar field perturbations around the real Choptuik solution does not add 
any unstable modes. The Choptuik solution is therefore an attractor of co-dimension one even for the complex scalar 
field. 

It would be interesting to run collapse calculations for a one-parameter family of matter models, such as the p = kp 
family |?J or the non-linear sigma model g] , where the Lyapunov exponents of the perturbations change continuously 
with the parameter. One might thus be able to find parameter regions where two equally strong attractors coexist, and 
new interesting phenomena would arise. It would also be extremely interesting if one could find a way of calculating 
or estimating the number of unstable modes of a given self-similar solution other than by constructing the solution 
and all its perturbations explicitly. 

Another point should be addressed in future theoretical work. As I mentioned in the introduction, the phase space 
picture of universality, although apparently correct in some aspects, is strictly speaking wrong because the same 
spacetime corresponds to many different trajectories in superspace, according to the way it is sliced. At the very 
least one should be able to derive a universal geometrical prescription fixing the lapse and shift from Cauchy data in 
superspace, such that the intuitive phase space flow is a realized as a Hamiltonian flow. Furthermore, this Hamiltonian 
flow should admit a geometric interpretation as a scaling, or "renormalisation group" , flow. The evolution in the time 
variable r (at constant £) does in fact go from one set of Cauchy data to the same, only changed in overall scale. 

Clearly the idea of approximate scale invariance is at the core of renormalisation group ideas and methods. Further- 
more, the calculation of 7 as given in section IV B is identical with the calculation of the critical exponent governing 
the diver gen ce of the correlation length given in any textbook on critical phenomena in statistical mechanics, for 
example [ f27| . As far as I can see, however, the "critical phenomena" analogy is not with critical phenomena in a 
system which is described by a partition function (such as a system in thermal equilibrium or a quantum field). There 
is neither a nonvanishing Hamiltonian, nor an inverse temperature /3, or quantum of action ft, such that one could 
construct a weight on phase space. 

The analogy seems to be rather with the application of renormalisation group methods to deterministic PDEs 
p8[ . This takes up ideas of Barenblatt's ^| of generalizing a self-similar ansatz such as f(r 2 /t) for a PDE such 
as a generalized diffusion equation to t a f(r^/t). Here, a is a non- integer "anomalous dimension" that cannot be 
determined by dimensional analysis. But the factor t a gives rise to terms proportional to a in the equation for /, and 
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this equation will admit regular solutions only for certain values of a: one has a nonlinear eigenvalue problem. 

In the critical collapse case we see this twice. In our ansatz for the perturbations we clearly have an explicit factor 
t x . But A has the same function in the background solution! There we expand all fields, which are periodic in r, in 
terms of e muJT f n ((), where u> = 2ir/A. In consequence we obtain terms proportional to u> in the equations for the 
expansion coefficients f n (()- The CSS case corresponds to these terms being absent. We may therefore consider the 
DSS ansatz as a generalization of the CSS ansatz parallel to the anomalous dimension in Barenblatt's examples, and 
consider iuj a complex anomalous dimension. These parallels certainly put the investigation of critical gravitational 
collapse into a wider context, and may yet give rise to new predictions or a simplified theory. 
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APPENDIX A: BACKGROUND BOUNDARY CONDITIONS 

£ = — oo and Q = are singular points of the equations. The boundary conditions are therefore implemented by 
expanding the field equations in powers of around £ = — oo and of £ around £ — 0. The resulting equations are 
given below. 

The regularity conditions at ( — ► — oo can be solved in closed form. We expand in powers of e 1 ", as 

o(C, t) = a {T) + ai(Y)e c + a 2 (r)e 2C + ... (Al) 

It turns out to be more natural to expand X and Y instead of X + and X_. As discussed in section ||, we impose 
do = 1 and go — 1. Expanding the field equations, we find that the two (periodic) functions Yi(rj) and £o(v) can be 
chosen freely. Their significance here is the following. £o parameterizes the class of spacetime coordinates we use, 
while a combination of the two contains free boundary data for the scalar field (f> f^]: 



-£oO)y- l( - r ) = 



<9(ln£) 

The other nonvanishing expansion coefficients up to order e 3 ^ in X and Y and order e 2l » in a and g are 



(A2) 



a 2 = -.92 = \y?, (A3) 
X 2 = \e^{Y>-{l + &)Y x ), (A4) 



3 

"1 



Y 3 = -^Y? + et° 



2 X^-(l + QX 2 



(A5) 



while the other coefficients vanish. These coefficients were calculated from ao = go = 1 and the evolution equations 
alone, but they also obey the constraint order by order. 

Regularity conditions at the boundary £ = give rise to an ODE system. We again make a power-law ansatz of 
the form 

o(C, t) = oo(t) + o x (t)C + a 2 (r)C 2 + ... (A6) 
Here we find that suitable fields to expand in are a, X + , A_ and, instead of g, the quantity 

D = (l + C )e^°g. (A7) 

Its evolution equation is easily seen to be 



4 This quantity was called Yq(t) in |16 



1G 



D x = D{2-a 2 ). 



(A8) 



The expansion coefficients can be caiculated recursively, by the solution of first-order ODEs, if we use X + o(t) and 
£o(t) as the free parameters. Our coordinate condition (24) is simply 



D = 1. 



We note that A_ drops out of the constraint (|T^) to leading order, which is 



a o = (l + fo)ao 



2 +a H^°-2 



(A9) 



(A10) 



Given A + o, this is a nonlinear, inhomogeneous, first-order ODE for ao. The one integration constant is fixed uniquely 
by the requirement that the solution be periodic. The regularity condition (|2^) gives us another ODE, this time 
linear, for X_ : 



We can now calculate algebraically, 



D x = 2 - at, 



1 



Xl 



A_ + (1 + £ )*+o - 0. 



«i = - fl o) + «o (^+o + Xlo)] » 



a nX+o 



X' 



+ -x +0 -x_o + (i + Q- i x' +0 



(All) 

(A12) 
(A13) 

(A14) 



For X-i we obtain again a linear ODE 



1 



X-x 



2aoaxX^o [ ^ + -^+0 ] + 2<2qA_oA + oA + i + X + x 



+ (2 - a 2 )X'_ = 0. 



(A15) 



To quadratic order, we have three algebraic expressions and one more linear ODE, which we do not give here. (We 
need the previous equations to determine X-x, which will be needed in equation (|B6| ) below.) We have used explicitly 
only the zeroth order of the constraint ( |l9| ) , but the first and second orders are satisfied identically as expected. 

Using these expansions we can calculate Cauchy data for a very large and a very small negative value of £, thus 
avoiding the vanishing of X± at ( = — oo and the breakdown of the Cauchy scheme at ( = 0. 



APPENDIX B: BOUNDARY CONDITIONS FOR THE PERTURBATIONS 

The perturbed boundary conditions at £ — > — oo in terms of the free perturbation 6iYx(r) are 



SiX 2 = ^e«° (SiY{ + (Aj — 1 — QSiYx) . 



S t Y 3 = -TYl^Yx + e io 



l -hX' 2 + (1a< - 1 - &)6iX 2 



At £ = the perturbed constraint (70) simplifies and no longer contains <5iAT_: 

1 



hao, T + AiJiao = (1 + £ ) 



1 + 1^(2X^-1) 



Sido + 2a^X +0 5iX +0 - -g L a^{X^ - Xi Q )5 igo 



(Bl) 
(B2) 
(B3) 

(B4) 



We can spec ify 5jX+g and Sigo freely and solve this equation for <5;a + o- The average value of the coefficient of Sia® 
in equation (B4) is \ + 1. When this vanishes, the equation has no solution (with periodic boundary conditions), 
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as discussed in the appendix on Fourier methods. This means that for A, = — 1, there are no perturbations that are 
regular on £ = 0. In the A plane this gives rise to a pole at A = — 1, see appendix (To calculate the average value 
of the coefficient, we note that the background constraint (fL9|), evaluated at £ = with the boundary condition 
reduces to 



(Inao). r = (l + Q 



(1 



4) + alX 2 +0 



As the left-hand side is the derivative of a periodic function, the right-hand side has vanishing average.) 
The vanishing of the numerator of the 8iX- equation is an ODE that can be solved for SiX-o, namely 



6iX' 



A,--(1 +;,',) (i(l-aS)-4*£ 



SiX_ + (1 + Q ao^-o (1 + 2Xl ) S iao 



+ ((1 + C^r 1 ^ - X-i) g^Sm + (2a 2 X +0 X_ + l) S t X +0 = 



(B5) 



(B6) 
equals 



This equation in turn has no solution if the coefficient of 5iX_o has vanishing average. This is the case if Ai 
the average of (1 + £q)(1 — a 2 ,), which numerically has value ~ —0.385. This gives rise to anoth er p ole. 

We do not need to expand the perturbations away from £ = 0, because in the discretisation (Dl) of the linearized 
equations we do not need to evaluate the £-derivatives of the perturbations at £ = 0. 



APPENDIX C: PSEUDO-FOURIER METHOD 



The discrete Fourier transform of the N complex numbers /„ is defined by 

N-l 

N 



1 N ~ X 



(CI) 



and its inverse is 



N-l 

fn = ^2 fkC 

k=0 



(C2) 



The motivation of this definition is of course that the f n are to represent the values of a smooth periodic (complex) 
function f(x) at N equidistant points over one period. The essential idea of pseudo-spectral (here: pseudo-Fourier) 
methods is to carry out algebraic operations pointwise on the /„, and integration and differentiation on the fk, 
switching from one to the other with a Fast Fourier Transform algorithm. A detailed description of pseudo-spectral 
methods can be found in ]30| ]. I give here only the technical information necessary to specify my numerical method. 

We shall need to define various operations on the Fourier components fk which are to represent operations on f{x). 
To do this in a consistent way, we have to start from a definition of f{x) in terms of the fk- We choose 



.V 



f(x) - f 



-1 

E 

k=l 



(fk( 



+ fN-kC' 



, Nir 

+ JK COS ( —X 



(C3) 



where A is the period of f(x). This expression obeys f(x n ) = fn for x n = nA/N, but this requirement alone does not 
uniquely define it. We can now derive how any operation on the fictitious function f(x) is represented as an operation 
on the fk, for example interpolation of f[x) to arbitrary values of x, differentiation or integration. Let us begin with 
differentiation. The corresponding operator D is given by 

D f) k =-^f>°> 0<k< N/2-1, 



N/2 



2ir(k — N) 



fk, N/2 + l<k<N-l. 



(C4) 
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We see that on differentiation the high-frequency cosine, which takes values of ±1 on the x n , is annihilated, because 
the high-frequency sine takes the value on the x n . For the purpose of integration we define a more general operator, 
namely I\ given by 



M 

(4- 



0, A = 0, 
1 



A - 2tt/c/A 

Ix f) »r,n = \2 _ t^AT/ A\2^/2 



/ fc , l<k< N/2-1, 



n/2 A 2 - (ttTV/A) 

(C5) 



Clearly, for A = and fo — 0, this is simply integration, with the integration constant fixed so that the integral has 
vanishing zero-frequency component. If f(x) has a nonvanishing average (zero-frequency part) fo, then its principal 
function has a part fox and is no longer periodic. In this case we enforce periodicity by setting fo equal to zero. 
When A ^ however, even for fo ^ 0, I\f solves the ODE 

(hf)' + Xhf = f, (C6) 

with the integration constant chosen such that I\f is itself periodic. With the help of this definition we can write the 
unique periodic solution / of 

/' + gf + h = (C7) 

in closed form, namely as 

/ = -e-^I- o ( e ^h) . (C8) 
This is of course only the standard use of an integration factor, written in Fourier components, and so that it obeys 



periodic boundary conditions. The expression diverges as go —* 0, and indeed the equation (07) has no solution with 
periodic boundary conditions for go = 0. All the ODEs we have to solve are of the form (C7) except the constraint 
(|l9|), which is nonlinear, of the form 

2a' = ga + ha 3 . (C9) 



Fortunately, it can be reduced to the linear form (C7) by the substitution / = a . 

We also need to define an operation that doubles N while representing the "same" function f(x), and its inverse. 
These operations will be required for "dealiasing" - doubling the fk before going to the /„, carrying out the necessary 
algebraic operations on the doubled /„ and going back to the /fc, then halving the fk and thus throwing away high 
frequency noise. With our definition of f(x), doubling means splitting the high frequency cosine. The doubling 
operation B from N to 2N components is 

Bf) =f k) 0<k< N/2-1, 



(Bf 



Bf) =0, N/2 + 1 < k < 3N/2 — 1, 
(Bf) k = f k - N , 2.N/2 + 1 < k < 2N - 1. 



(C10) 

The halving operation S from N to N/2 components is 
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k 



(Bf) 



fk, < k < N/4 - 1, 

/at/4 + f3N/i, 



N/4 

Bf) ; = /fe+Ar/2, A/4 + 1 < k < N/2 - 1. 

(Cll) 



I stress once more that the definitions of D, I\, B and S follow uniquely from the definition (C3). With our choice, 
D and Iq are not the inverse of one another because of the treatment of the high-frequency cosine. SB is the identity 
operator, while BS is a smoothing operator. Our choice of definition is motivated by the requirement that none of 
the operations mix the real and imaginary part, that is Df , I\f, Bf and Sf are all real if / is real. I require this 
because I typically encode two independent real functions as the real and imaginary parts of one complex function. 

In a previous code |jl6), I did not use one consistent definition of f(x) in the Fourier discretisation. More seriously, 
the discretisation mixed real and imaginary parts in some places even where the real and imaginary part represent 
two unrelated real functions. I have therefore replaced it by the one described above. This results in a small change 
in the numerical results for the same equations discretized on the same grid, in particular a slightly changed value of 
A. 



APPENDIX D: NUMERICAL METHODS FOR THE CRITICAL SOLUTION 



I use Fourier transformation to transform 1 + 1 hyperbolic PDEs into an ODE system, and ODE boundary conditions 
into algebraic boundary conditions. I use the assumption that X± are odd while a, g and £o are even in r to encode 
the variables X±, g and £o as a single complex one, as £o + ig + X + + iX-. The discrete Fourier transform of this 
function constitutes my independent variables. Let us assume that the information contained in them is to correspond 
to a sampling of X±, g and £o at N points each. Then the total number of complex variables is also N, or 2N real 
variables. Typically I use 2N = 128. Out of these 2N variables, I set the one combination corresponding to the 
fundamental sine mode of £ equal to zero, storing A in its place. As mentioned in the main text, this fixes the 
translation invariance in r of the problem and balances the degrees of freedom against the constraints. In the analytic 
continuation problem the complex variable function is T + iG + X + + iX-. For the linearized perturbations it is 
iSig + 6iX + + iSiX_, plus the constant A^, corresponding to 3 A/2 + 1 real variables. The linear perturbations are 
normalized by an additional condition at the right boundary, namely that the root-mean-square of X + q be 1. 

£o and A are formally part of the dependent variables of the background problem, although with vanishing £- 
derivative. In the perturbation problem they are supplied as external fixed parameters, together with the background 
solution. In each case, a or Sid is reconstructed at each step from the other three fields by the solution of the constraint. 
This solution can be given in closed form in Fourier components, as explained in appendix Previously, we have 
doubled the number of Fourier components, transformed back to real space, and separated the real and imaginary and 
the even and odd parts to obtain our four fields sampled at 2N points each. At the same time we have calculated the 
necessary r-derivatives. Now we calculate the necessary algebraic expressions pointwise in the 2N sampling points, 
put the result back together in complex form, transform back, and halve the number of Fourier components. 

The ODE system df /d( — F(f, £) is discretized in a standard way as 



f 



n+1 



f n — {Cn+1 - Cn) F 



fn+1 + fn Cn+1 + Cri 



(Dl) 



To solve this system of algebraic equations to machine precision together with algebraic boundary conditions on two 
sides, I use a standard relaxation algorithm |nj. If I have to integrate these equations forward (that is, with boundary 
conditions ("initial data") only on one side), I use a second-order Runge-Kutta step as a first guess for f n +i, and then 
solve the discretized equations to machine precision using Newton's method. On one hand this is necessary because 
an explicit algorithm is unstable for these equations even for very small time steps. On the other I want the equations 
to be discretized in the same way as in the relaxation algorithm when I continue the solution from £ = to large C- 
Schematically, F(f, — N(f, ()/D(f, Q. As discussed in the main text, the SSH announces itself by D = 0. As a 
boundary condition we impose that D = at £ = for al l r. We then impose N = as a second boundary condition 
at C = 0. C — is itself the last grid point. We see from (Dl) that F is never evaluated at D — 0. There is therefore 
no need to treat the singular point £ = in a special way. If we want to shoot away from Q — 0, however, for example 
in order to prime the relaxation algorithm, we need to expand the field equations in powers of (. This is done in 
appendix |a| 
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The other set of boundary conditions are imposed at £ = — oo. We do this by expanding in powers of exp£ (see 
appendix . 



APPENDIX E: NUMERICAL TREATMENT OF THE PERTURBATIONS 

I have studied the perturbations with two largely independent codes. One code assumes that Xi is real and that 
the perturbations have the same symmetry (X± odd and a and g even) as the background. A second code allows 
for complex Xi and perturbations, and the perturbations either having the same symmetry as the background, or 
the opposite symmetry. The two kinds of perturbations decouple because of the background symmetry and can be 
considered separately. 

I used the following discretisation of the linearized version d5f /d( = dF(f, ()/df ■ 5f of the ODE system df /dC, = 
F(f,()- 

( \idF\~ 1 ( hdF\ r , s 



2 df J \ 2 df 

The matrix dF(f, / df is extracted in Fourier components, by varying one Fourier component of Sf at a time in the 
linearized derivatives dF(f, Q/df -8f . Its coefficients are evaluated on the background taken at the midpoint between 
/„ and f n +i, and h = Cn+i — Cn- This scheme is explicitly linear, implicit (and stable), and second-order accurate. 

For finding the eigenvalues A^, both codes generalize ideas already used in similar calculations |To[ ], using the 
linearity of the equations. All possible linear perturbations compatible with the boundary conditions at either C = 
or £ = — oo are evolved together to a midpoint, say £ = — 1. One obtains a complex square matrix d (mismatch at 
the midpoint) / d(free parameters at the endpoints). As a test, I have varied the real and imaginary part of the free 
parameters independently. The Cauchy-Riemann equations (with respect to the free parameters) are obeyed without 
any numerical error. 

This matrix, say A, depends on A^. If for a given Xi a linear perturbation exists which obeys all boundary conditions 
at both endpoints, det A = 0, and this perturbation is the eigenvector with eigenvalue zero. We therefore have to search 
for zeros of det A(Xi) in the complex plane. Because A* does not appear in the equations, det A(Xi) is holomorphic. 
Complex conjugation of A corresponds to a certain interchange of both rows and columns (the interchange of positive 
and negative frequencies in r), and therefore det A(A*) = det A(Xi)* . I have checked both symmetries numerically. 
I find a relative numerical error in the Schwarz reflection principle of 1CP 10 and in the Cauchy-Riemann equations 
(with respect to A;) of 10~ 9 . 

It is sufficient to consider the strip < 5Ai < A/-7T. Because det A(Xi) is holomorphic, we can use the well-known 
contour integral J A'/AdXi — 2iti(N z — N p ) to count the number of zeroes minus the number of poles within the 



contour. Apart from the perturbations discussed in section IV, I find a simple pole at A,; = —1, for modes with the 
same symmetry as the background, and a simple pole at A ~ —0.385 for modes with the opposite symmetry. Their 
origin is discussed in appendix |B| 

Once I have found a zero of det A(A^), I determine the zero eigenvector, and hence the corresponding perturbation 
field. (Numerically I find the eigenvector by singular value decomposition.) To refine this result, and to provide 
an independent check, I use it as the starting point of a relaxation algorithm. In this algorithm, Xi is one of the 
dependent variables, and its presence is balanced by a boundary condition fixing the norm of the linear perturbation. 
I have implemented this algorithm only for real Xi, as I am mainly interested in the unstable mode, which is real. 

APPENDIX F: NUMERICAL ERROR ESTIMATES 



For an ODE discretisation to second order such as (Dl) one would expect the solution to converge globally to 



second order in the step size. I had incorrectly claimed to observe this in jig] . The real situation is more complicated: 
the root-mean-squared or maximum measure of convergence show convergence that is quadratic on the average over 
a wide range of step sizes. By quadratic on the average I mean that the log-log plot of step size versus the difference 
of two numerical solutions with different step sizes is a somewhat wiggly line, but with an average slope of two, as 
illustrated in Fig. |3[ This result is puzzling, because one would expect a much tighter result to hold, namely that 
the numerical solution of the equations at a given step size h is equal to the continuum solution plus a discretisation 
error which itself is a power series in even powers of h only (the Richardson form). This clearly does not hold here. 
I believe that the irregularity is due either to the presence of boundary conditions, or in particular to the singular 
nature of the boundaries. 
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The fact that the numerical solution really obeys (Dl) to machine precision is easily established, independently 
of the (complicated) relaxation algorithm that found the solution in the first place. Furthermore the convergence 
behavior should be independent of what is inside the routines that supply the algebraic boundary conditions and the 
derivative F(f, £). Therefore it appears unlikely that the cause of the irregular convergence behavior is a programming 
error, even if such an error was present. 

When a numerical error has the Richardson form, one obtains a sharp error estimate, and which moreover is based 
on a tested theory of the numerical solution process. Here I cannot use this theory. Nevertheless, there clearly is 
something like quadratic convergence with decreasing step size, and so I feel justified to use the difference between the 
numerical solution with N grid-points and the numerical solution with 2N grid-points as a measure of the numerical 
error of the solution with N gridpoints. Adding a safety factor of 4, I take this value as the error of the 2N solution, 
and the 27V solution itself as the best value. 

Fig. |3| plots maximum and root-mean-square measures of convergence of the entire solution. As it happens, the 
maximum measure is dominated by the error in A, and so doubles as an error plot for A. As described above, I 
obtain the value and error bar A = 3.4453 ± 0.0005. 

Table || gives the values of the parameters A, Ai and A2 for various N. A2 is the gauge mode corresponding to 
a translation of the background solution, and we should find A2 = 0. Its absolute numerical value gives a bound 
on the error of A of ±0.009, while from the convergence of Ai I would have estimated a much smaller error of 
±0.0002. I opt for the larger error estimate, and obtain Ai = —2.674 ± 0.009, corresponding to a critical exponent 
7 = -1/Ax = 0.374 ±0.001. 



APPENDIX G: CURVATURE AT THE FUTURE SSH 

The Riemann tensor in spherical symmetry in general coordinates has been calculated in p2]| . Substituting the 
form (Q) of the metric, it is straightforward to verify that the only appearance of any second derivatives of a or a 
with respect to r and t in any component of the Riemann tensor is in the combination 

±r(r&) -r(r^) . (Gl) 



a \ a / ,t \ a 
We now use the Einstein equations 

= \a> (Xl A 2 ) , = i(a 2 - 1) + i« 2 (x\ - Xl) , (G2) 

to transform the terms in round brackets into algebraic expressions, and we transform the outer derivatives from 
coordinates (r,t) to coordinates (w, p) via 

Q Q Q Q ^ Q 

^ dt dw' Or dp W dw 

The only derivative in the resulting expression that potentially does not exist is A + ul . But it arises as 

{G + Tw){a 2 X 2 + ), w , (G4) 

and while X+ >w blows up as (w + l) 6 " 1 as w — ► —1, the coefficient vanishes as (w + 1), because of the coordinate 
condition T = G at w = — 1. We conclude that the curvature remains finite at w = — 1 although not all second 
derivatives of the metric are finite there. 
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TABLE I. A comparison between the Kerr-Newman and critical collapse solutions. 



Kerr-Newman solutions 



Critical collapse solutions 



length scale M 
stationary 

horizon 
==> vacuum 
quasi-normal modes e' w *'' M 
attractor 
"no hair" 



scale-invariant 
self-similar 
naked singularity 
=> matter-dependent 
perturbations t Xi 
attractor of co-dimension one 
=>■ universality 



TABLE II. Convergence of A and A with step size in Delta is the echoing period. Ai is the Lyapunov exponent of the 
one growing mode. Its negative inverse is the critical exponent 7. A2 is the exponent of the translation gauge mode. It must 
be zero and serves as a check on the numerical error. Note that the numerical grid (and number of steps) is the same for the 
background as for the perturbations in each case. 



Number of steps A Ai A2 

51 3.4321725669119 -2.6858281957399 -3.8648655101422D-02 

101 3.4513015765429 -2.6700545097384 6.3598135051540D-02 

201 3.4431664827331 -2.6741157762787 -2.4122047436405D-02 

401 3.4446384162424 -2.6748914760819 -2.3339465009179D-02 

801 3.4458770431665 -2.6738878803912 2.9009984141595D-02 

1601 3.4453484479734 -2.6740958070987 -9.4268323369794D-03 
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FIG. 1. The phase space picture for discrete self-similarity. The plane represents the critical surface. (In reality this is a 
hypersurface of co-dimension one in an infinite-dimensional space.) The circle (fat unbroken line) is the limit cycle representing 
the critical solution. The thin unbroken curves are spacetimes attracted to it. The dashed curves are spacetimes repelled from 
it. There are two families of such curves, labeled by one periodic parameter, one forming a black hole, the other dispersing to 
infinity. Only one member of each family is shown. 
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FIG. 2. Schematic diagram of the coordinate systems I use. The horizontal and vertical axes are r and —t. Horizontal 
accumulating lines are r = 0, r = A, r = 2A, etc., from bottom up. Vertical accumulating lines are p = 0, p — A, p — 2A, etc., 
from right to left, a) The curvature singularity (r = 0, t = 0). b) The regular center r = 0. c) The past self-similarity horizon 
C = 0. d) The matching line £ = £o, w = 1. e) The matching line w = 0, t = 0. f) The future self-similarity horizon w — — 1. 
I use four separate coordinate patches: 1) The nonlinear eigenvalue problem for A. 2) Continuation to Co- 3) Continuation to 
t = 0. 4) Continuation to the future SSH. 
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FIG. 3. Convergence of the solution from r = to the past SSH. Horizontal axis gives number of grid points (26, 51, 
801). Vertical axis gives difference of numerical solution from reference solution with 1601 points. Upper line is maximum error 
over all components and grid points, lower line is root-mean-square error. Upper line coincides with the error in A, as this is 
the component which has the largest error. Note that the diagonal of the box has slope 2, the expected convergence behavior. 
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